#Pop 1 vs Pop 2 Fst
#try to find SNPs associated with high/low FEC
setwd("/Paterson/Datafiles/grouse/grouse_analysis")
load('wksp_2011_01_11.rdata')

grouse.df4 <- grouse.df3[,-grep('pop.3',names(grouse.df3))]
grouse.df4 <- grouse.df4[rowSums(grouse.df4[,c('pop.1.freq','pop.2.freq')])>0,]

grouse.df4$p.hat <- apply(X=grouse.df4,1,function(X){
	#X <- as.vector(X)
	#class(X)
	#X[grep('freq',names(X))]
	sum(as.numeric(X[grep('freq',names(X))])*as.numeric(X[grep('depth',names(X))]))/sum(as.numeric(X[grep('depth',names(X))]))
	})

ifst <- seq(from=0.001,to=0.99,length.out=100)
itheta <- 1/ifst -1
lik.theta.fit2 <- matrix(0,nrow=nrow(grouse.df4),ncol=100,dimnames=list(grouse.df4$id,paste("theta",itheta,sep="_")))


for(i in 1:100){
	tst <- numeric(nrow(grouse.df4))
		for(j in 1:length(tst)){
		lik.theta.fit2[j,i] <- LikBetaBinomialPops(depth=as.numeric(grouse.df4[j,grep('depth',names(grouse.df4))]),
		count=as.numeric(grouse.df4[j,grep('count',names(grouse.df4))]),
		p.hat=as.numeric(grouse.df4[j,'p.hat']),theta=itheta[i])
		}
	#lik.theta.fit[i] <- sum(tst)
	}
null.which2 <- which(colSums(lik.theta.fit2,na.rm=TRUE)==max(colSums(lik.theta.fit2,na.rm=TRUE))) #4, Fst = 0.03

BBlikTest.df2 <- data.frame(null.lik = lik.theta.fit2[,null.which2],max.lik = rep(0,nrow(grouse.df4)))
BBlikTest.df2$max.lik <- apply(lik.theta.fit2,1,max)
BBlikTest.df2$chi.sq <- 2*(BBlikTest.df2$max.lik-BBlikTest.df2$null.lik)
BBlikTest.df2$p.val <- 1-pchisq(BBlikTest.df2$chi.sq,1)

library(multtest)
tmp <-  mt.rawp2adjp(BBlikTest.df2$p.val,'Bonferroni')
BBlikTest.df2$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']

tmp <-  mt.rawp2adjp(BBlikTest.df2$p.val,'Hochberg')
BBlikTest.df2$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']

#repeat with a dodgy method to pool snps at each locus
BBlikTest.pooled2 <- aggregate(lik.theta.fit2,by=list(grouse.df4$cDNA),sum)
BBlikTest.pooled2.df <- data.frame(cDNA=BBlikTest.pooled2[,1],null.lik = BBlikTest.pooled2[,null.which2+1],
	max.lik = rep(0,nrow(BBlikTest.pooled2)))
BBlikTest.pooled2.df$cDNA <- as.character(BBlikTest.pooled2.df$cDNA)	
BBlikTest.pooled2.df$max.lik <- apply(BBlikTest.pooled2[,-1],1,max)
BBlikTest.pooled2.df$chi.sq <- 2*(BBlikTest.pooled2.df$max.lik-BBlikTest.pooled2.df$null.lik)
BBlikTest.pooled2.df$p.val <- 1-pchisq(BBlikTest.pooled2.df$chi.sq,1)

tmp <-  mt.rawp2adjp(BBlikTest.pooled2.df$p.val,'Bonferroni')
BBlikTest.pooled2.df$bonferroni <- tmp$adjp[order(tmp$index),'Bonferroni']
tmp <-  mt.rawp2adjp(BBlikTest.pooled2.df$p.val,'Hochberg')
BBlikTest.pooled2.df$hochberg <- tmp$adjp[order(tmp$index),'Hochberg']
BBlikTest.pooled2.df$fst <- 0

BBlikTest.pooled2.df$fst <- apply(BBlikTest.pooled2[,-1],1,function(X){
	fst.seq <- seq(from=0.001,to=0.99,length.out=100)
	fst.seq[which.max(X)]
	})

BBlikTestSigGenes2 <- grouse.df4[grouse.df4$cDNA %in% BBlikTest.pooled2.df[BBlikTest.pooled2.df$hochberg<0.05,]$cDNA,]

save.image('wksp_2011_03_08.rdata')